num_cores = str2num(getenv('SLURM_CPUS_PER_TASK'));
c = parcluster;
poolobj = parpool(c, num_cores);

clear all
clc

cd('../cvx_2.2.2')
cvx_setup

warning('off','MATLAB:nargchk:deprecated')
warning('off','MATLAB:illConditionedMatrix')

cd('../experiment')

data = readtable('data.csv');

J = 3;    % # of goods

data_wide.choice = data.choice;
data_wide.x = [data.x1 data.x2 data.x3];
data_wide.z = [data.z1 data.z2 data.z3];

data_wide = struct2table(data_wide);

[data_trans] = transform_data_01_exp(data_wide,J);

m_regr = [3 4 3*ones(1,J-2) 4 4 3*ones(1,J-2)];

[combos_all_unique] = combinations_10(m_regr,J);

% constraints

[Aineq,bineq] = constraint_mat_theta_2(combos_all_unique,J);

Aeq = [];
beq = [];

ntheta = size(Aineq,2);

N = size(data_wide,1);

Nbs = 1000;

ratio_trimmed_mean_bs = NaN(Nbs,1);
mean_ratio_bs = NaN(Nbs,1);
trimmed_mean_ratio_bs = NaN(Nbs,1);
theta_NPD_bs = NaN(Nbs,ntheta);

for bs=1:Nbs
    
    % draw random numbers for bootstrap sample (we ran it in multiple batches, which is why the seed for the random numbers has the expression below)
    if bs/Nbs <=1/5 && bs/Nbs>0
        temp = bs;
        rng(temp)
    end
    for i=2:5
    if bs/Nbs <=i/5 && bs/Nbs>(i-1)/5
        temp = (bs-(i-1)*200)*i*1000;
        rng(temp)
    end
    end
    index = randsample(N,N,true);
    
    data_trans_bs = cell(J,1);
    for j=1:J
        data_trans_bs{j} = data_trans{j}(index,:);       
    end
    
    [theta_NPD,fval_NPD,flag_NPD] = estimation_NPD_cvx(data_trans_bs,J,N,m_regr,Aineq,Aeq,bineq,beq,...
        combos_all_unique,ntheta);
    
    theta_NPD_bs(bs,:) = theta_NPD;
    
    [ratio_trimmed_mean_bs(bs),mean_ratio_bs(bs),trimmed_mean_ratio_bs(bs)] = script_ratio_rand_exp_fun_average(theta_NPD,m_regr,data_wide,J,ntheta,combos_all_unique);
    ratio_trimmed_mean_bs(bs)
    save('bootstrap.mat')
    bs
end

std(ratio_trimmed_mean_bs)

clear data* A* b*  der_den der_num
    
save('bootstrap.mat')
